论文推荐 | 黄炎, 王庆宾, 冯进凯, 等：局部地形改正快速计算的GPU并行的棱柱法

《测绘学报》

构建与学术的桥梁        拉近与权威的距离

复制链接，关注《测绘学报》抖音！

【测绘学报的个人主页】长按复制此条消息，长按复制打开抖音查看TA的更多作品##7NsBSynuc88##[抖音口令]

本文内容来源于《测绘学报》2020年第11期，审图号GS（2020）5931号。

局部地形改正快速计算的GPU并行的棱柱法

黄炎¹, 王庆宾¹, 冯进凯¹,邢志斌², 范雕¹, 谭勖立¹,吕明昊¹                            

1.信息工程大学, 河南 郑州 450001;
2.航天工程大学, 北京 102200

基金项目：国家自然科学基金（41574020）；信息工程大学校自立课题（2105070232）

摘要：在重力归算中，局部地形改正在重力勘探、地壳结构分析和大地水准面计算等领域有着重要意义，但严格棱柱体积分公式计算效率低，而快速计算公式则会降低计算精度。本文利用CUDA并行编程平台，提出一种地形格网重新编码和严格棱柱体积分八分量拆解方法，实现了基于CPU+GPU异构并行技术的严格棱柱体积分计算地形改正快速并行算法，克服了GPU各个线程计算任务分配和线程计算超载问题，解决了局部地形改正的高分辨率、高精度严密公式的快速计算难题。通过试验，在显卡型号为Tesla V100的计算机上进行4°×6°范围，积分半径40'和分辨率1'的局部地形改正计算仅需1.5 s；分辨率10″的局部地形改正计算仅需14.6 min；进行分辨率3″的地形改正计算耗时45.7 h，而传统串行算法则难以完成计算。在保证微伽级以上计算精度的条件下，计算加速比最高达到850倍以上，有效缩短了计算耗时，提高了计算效率。本文还依据上述并行算法对全国范围地形改正量进行计算。结果表明，我国地形改正量普遍低于80 mGal（1 Gal=10^-2 m/s²），平均值1.83 mGal，最大值达到196 mGal。

关键词：地球重力场    局部地形改正    GPU    CUDA    严格棱柱体积分    加速比  

引文格式：黄炎, 王庆宾, 冯进凯, 等. 局部地形改正快速计算的GPU并行的棱柱法. 测绘学报，2020，49(11)：1430-1437. DOI: 10.11947/j.AGCS.2020.20190403.

阅读全文：http://xb.sinomaps.com/article/2020/1001-1595/2020-11-1430.htm

全文概述

重力归算是重力场资料解释与应用的重要工作^[1-3]。为计算方便，通常将归算分为3部分：空间改正、层间改正和局部地形改正。其中局部地形改正通常基于棱柱法计算，其计算耗时随着地形分辨率和积分半径的增加而呈指数增加，计算效率低。例如，使用30″×30″的地形数据，积分半径取2°，在CD-4360计算机上，利用严格棱柱体积分计算一幅4°×6°范围区域局部地形改正，需要耗时384 h(文献[4])。为提高计算效率，通常将计算区域分为中央区和近区。中央区采用棱柱法，近区则采用线性卷积积分法进行逼近。尽管计算效率得到了提高，但精度有所损失。国内外学者在保证一定精度的前提下，利用快速傅里叶变换(FFT)进行地形改正计算，对4°×6°范围、30″分辨率，积分半径为2°的区域进行局部地形改正计算，效率提高近100倍。但在地形起伏剧烈区域精度有所损失，在我国高山地区误差最高达到20 mGal^[5-9]。因此，如何在不损失精度的前提下提高计算效率成为了研究重点。随着计算机技术与大数据运算的发展，大规模数据运算不仅仅局限于专业且昂贵的并行计算机，GPU的快速发展使得大规模并行运算成本急速降低^[10-15]。GPU的线程数量和计算能力是CPU数十倍乃至上百倍，充分利用GPU的计算能力成为并行计算的关键。本文利用CUDA运算平台，基于CPU+GPU(显卡)的异构并行模式，构建了棱柱法计算局部地形改正的异构并行算法公式，充分利用CPU的调控能力和GPU的计算能力^[10]，解决了高精度、高效率局部地形改正计算问题，并对其计算精度和效率进行分析。

1  地形改正计算方法

1.1  严格棱柱体积分法

以待求点P在高程起算面的投影O为原点建立坐标系，x轴指向东，y轴指向北，P_i为待积分格网点，如图 1所示。图中，P表示待求点，(x_p, y_p, h_p)表示P在格网中的水平坐标和其高程；P_i表示待求点周围积分半径内的待积分格网点，(x_i, y_i, h_i)表示P_i在格网中的水平坐标和其高程。

图 1 地形格网棱柱 Fig. 1     Topographic grid prism diagram      

图选项

根据局部地形改正的定义，在图 1所示的坐标系下，积分区域对单个地形网格单元产生的局部地形改正值算法公式为

 (1)

式中, δg为积分区域对单个地形网格单元产生的局部地形改正值；σ表示积分区域；G为万有引力常量，G=6.754×10^-11 N·m²/kg²；ρ为地壳密度均值，取地球表面岩石密度近似表示地壳密度，ρ=2.67 g/cm³。

为了更好地进行编程计算，将上述公式有关物理量根据地形格网数据要求，以格网形式进行离散化，可得δg的离散化积分表达式如式(2)

 (2)

式中，x_p、y_p分别表示待求点在格网中的横纵坐标；x_i、y_i分别表示待求点周围积分半径待积分点在格网中的横纵坐标；h_p、h_i分别表示待求点与周围点的高程值；Δx、Δy为平面格网近似步长，Δx=RΔB、Δy=Rcos B_mΔL，R为地球平均半径，R=6 378 137 m，B_m为球面区域平均纬度，ΔB、ΔL为球面格网分辨率；s_L、s_B分别表示经度和纬度方向的积分半径点数。通过上述格网化离散积分公式，便可利用严格棱柱体积分编程计算局部地形改正δg。

1.2  基于CUDA的棱柱法积分并行算法

尽管将计算区域分为中央区和近区，中央区采用棱柱法，而近区采用线性卷积积分法可有效提高计算效率，但计算精度难以保证。尤其是在地形起伏较大区域，线性卷积积分相较于严格棱柱体积分误差会增大，高山地区误差最高可达20 mGal^[7-9]。所以，一直以来，局部地形改正计算的精度和效率难以同时满足。为有效提高棱柱法局部地形改正计算效率和精度，本文提出利用CPU+GPU的异构并行模式对上述局部地形改正算法进行并行改进。首先，依据GPU计算特性并结合严格棱柱体积分的计算过程对其进行并行策略设计；然后在可行的并行策略的支撑下，依据GPU线程计算能力，对原有严格棱柱体积分算法进行改化，使之能够更好适用于GPU计算逻辑。

1.2.1  GPU并行策略

利用严格棱柱体积分计算局部地形改正，就要计算地形格网中每一个点周围积分半径内的格网点对其的影响，如图 2所示。其中，n、m表示计算区域格网行数和列数，r表示积分半径所包含的地形格网个数(即积分半径/地形格网分辨率)。

图 2 局部地形改正索引 Fig. 2     Indicator sketch of local topographic correction      

图选项

而利用CUDA对上述过程进行并行化，需要着重解决两个问题，一是在开辟线程格网时，尽量使用较低维度的线程格网，可以有效减少线程调度过程的耗时^[16-19]；二是如何在低维格网中对对应的高维格网点进行索引。

(1) 格网降维。首先，确定并行计算单元(即GPU线程循环计算最小单元)。对于上述局部地形改正计算，可以以计算一个格网点地形改正量为计算单元，不同点位的计算分配给GPU的不同线程。在确定计算单元后，直接依据地形格网对计算单元进行分配需要开辟二维线程格网，这会导致GPU索引线程耗时增加^[20]。为有效减少线程索引耗时，减少GPU负荷，本文对计算区域格网进行重新排列，如图 3所示，将图 2计算区域格网以行为单位、首尾相接依次排列至一维数组中，数组长度为n×m。

图 3 局部地形改正GPU索引 Fig. 3     GPU index map for local topographic correction      

图选项

通过上述排列，原二维地形格网被转化为一维地形数据数组，长度为n×m，用tid表示待求点在一维数组中的位置编码。这样可以有效利用GPU的计算能力，且线程格网也仅需要最易检索的一维数组。

(2) 低维格网索引设计。依据上述格网降维方案，二维格网被有序的排列至一维格网(数组)中，因此对于待求点的索引较为简单，但待求点周围积分半径内的格网点则难以直接索引。为解决索引问题，本文引入row_ID数组，用于存储每个格网点的行号n_i，在进行并行计算时，仅需要将row_ID数组传入GPU端，即可快速得到待求点行号n_i=row_ID[tid]，通过行号便可计算列号m_i=tid-n×n_i，从而完成对待求点周围格网点的索引。通过上述并行策略，就可以根据一维线程块和线程，对需要计算局部地形改正的格网进行索引，利用CUDA的内建一维索引变量blockDim.x、blockIdx.x和threadIdx.x完成对地形格网的遍历。其中，blockDim.x表示每一个线程块所包含的线程个数；blockIdx.x表示当前调用线程所在线程块的编码；threadIdx.x表示当前调用线程在其所在线程块中的位置编号。通过这3个内建变量的组合，就可以完成对GPU线程的精确调度。具体索引如式(3)所示

 (3)

式中，tid表示待求地形格网点编号，范围为0~(n×m-1)。

1.2.2  严格棱柱体积分拆解

通过1.2.1节的并行策略，每个GPU线程平均分配若干个待求点进行计算，可以有效利用GPU的计算能力，对于低分辨率地形格网数据，可以直接采用式(2)进行计算。但是随着地形格网分辨率与积分半径的提高，每个GPU线程的计算任务与计算耗时呈指数增长，单个计算单元的计算量显著上升，会造成GPU线程计算超负荷而导致计算崩溃。经过分析，造成计算崩溃的主要原因是对f(x, y, z)的积分运算严重影响了GPU线程的计算速度，反复调用函数会增加GPU的负荷，因此还要对传统棱柱法公式进行改进。本文考虑将三重积分进行展开，并进行进一步拆解。

本文提出局部地形改正计算八分量法，即根据严格棱柱体积分过程将局部地形改正计算式(2)平均拆解为8个部分(三重积分，每次产生两个做差的式子，最终产生2³个式子)。如式(4)所示

 (4)

式中，i₁=x_i-x_p+0.5;i₂=x_i-x_p-0.5;j₁=y_i-y_p+0.5;j₂=y_i-y_p-0.5;Δh=h_i-h_p; v₁、v₂、v₃、v₄、v₅、v₆、v₇、v₈分别表示拆分后的8个分量。

利用上述棱柱体积分八分量计算式进行并行计算后，将不同格网点的各个分量按照图 3所示排列方式分别存储至GPU端的8个长度为n×m的数组中，将这8个数组分别表示为v₁[tid]、v₂[tid]、v₃[tid]、v₄[tid]、v₅[tid]、v₆[tid]、v₇[tid]、v₈[tid]。为减少数据传输产生的多余耗时，每一个分量计算完成后不释放GPU端显存，待所有分量均计算完成，无须将其传导至CPU端即可在GPU端对其进行并行向量求和，如式(5)所示

 (5)

式中，δg[tid]是一个长度为n×m的数组，用于存储各个格网点的地形改正量。利用上述式(4)、式(5)即可完成高分辨率局部地形改正的计算。

2  试验验证与分析

为验证本文方法的有效性，采用下述硬件、软件试验平台对该算法进行试验验证。

试验环境：高性能台式机，处理器为Intel® Xeon(R) Gold 6130 @ 2.10 GHz；GPU为Tesla V100；内存64.00 GB。软件环境为64位Windows Sever 2012操作系统，VS2013+CUDA10.1平台。

2.1  并行算法精度分析

为验证本文所提并行算法精度，选取4°×6°范围(地形高程格网范围为6°×8°。该区域最高高程3061 m、最低高程1 m、区域平均高程339.78 m；中间4°×6°区域最高高程2146 m、最低高程14 m、平局高程407.04 m)，区域高程具体如图 4所示。

图 4 区域地形 Fig. 4     Regional topographic map      

图选项

在CPU端，利用传统严格棱柱体积分串行算法计算图 4中间白框内4°×6°范围，分辨率为1′，积分半径为40′的局部地形改正，并以此结果作为“真值”。再分别利用CPU设备，近区域使用线性卷积分公式(中央区积分半径1′范围内使用严格棱柱体积分)计算该区域同分辨率、同积分半径的局部地形改正值；利用GPU设备，基于上述八分量并行算法公式，计算该区域同分辨率、同积分半径的局部地形改正值。将其结果分别与“真值”进行对比，对比结果见表 1。

表 1 不同计算方法精度分析Tab. 1 Accuracy analysis of various methodsmGal

计算方法	最大差值	最小差值	平均差值	RMS
严格棱柱积分八分量法	10-6	-10-6	< 10-12	< 10-12
近区使用线性卷积分法	4.77	-0.13	0.11	0.32

表选项

近区使用线性卷积分法与传统严格棱柱体积分串行算法的具体差值分布如图 5所示。

图 5 严格棱柱体积分与近区使用线性卷积积分差值 Fig. 5     Integral difference diagram of strict prism method and linear winding machine      

图选项

由表 1、图 5可以看出，八分量并行算法相较于传统严格棱柱体积分最大误差绝对值仅为10^-6 mGal，平均误差和均方根误差均小于10^-12 mGal，充分证明了本文所提八分量并行算法的精度可靠性。而在近区使用线性卷积分公式计算时，计算误差会随地形起伏的幅度变化，地形起伏剧烈地区误差较大、地形平缓地区误差较小，绝大多数点位误差保持在2 mGal以内，最小差值-0.13 mGal，个别点位差值较大，地形起伏剧烈区域误差最高达到4.77 mGal；平均差值0.11 mGal；RMS值0.32 mGal。尽管使用Toeplitz矩阵可以有效提高线性卷积分的计算效率^[21-23]，但精度损失依旧无法弥补。

2.2  不同分辨率地形格网计算结果分析

本节对比分析了不同分辨率地形格网数据计算局部地形改正的耗时，计算区域如图 4所示4°×6°范围。分别对分辨率为1′、30″、10″、3″，积分半径为40′的局部地形改正进行计算。其中，近区使用线性卷积积分进行计算时，中央区积分半径1′范围内使用严格棱柱体积分进行计算。结果见表 2。

由表 2可以看出，随着分辨率的增加，棱柱法计算局部地形改正的耗时呈指数增加，高分辨率地形改正计算难以在短时间完成，本文所提八分量并行算法可以有效提高棱柱法积分计算局部地形改正的计算效率。计算不同分辨率局部地形改正的加速比最高到达850倍以上，完成一个4°×6°范围1′分辨率地形改正计算仅需1.5 s；完成一个4°×6°范围10″分辨率地形改正计算仅需15 min，而传统棱柱法积分则需要一周以上时间。完成一个4°×6°范围3″分辨率地形改正计算用时45.7 h，而传统串行算法则难以完成。

表 2 各种方法不同分辨率计算耗时Tab. 2 Computing time-consuming of various methods with different resolutions

分辨率	待求点数	积分半径内待积分点数	严格棱柱体积分串行计算耗时	近区使用线性卷积积分计算耗时	严格棱柱积分八分量法并行计算耗时
1′	86 400	6560	9.682 min	62.518 s	1.453 s
30″	345 600	25 920	90.254 min	16.370 min	16.811 s
10″	3 110 400	231 360	208.266 h	21.940 h	14.633 min
3″	34 560 000	2 563 200	—	—	45.676 h

表选项

2.3  不同积分半径计算结果分析

本节对不同积分半径范围局部地形改正计算进行统计与分析，使用地形格网分辨率为10″。由于积分半径一般选取30~70 km范围^[24-25]，因此本节分别对积分半径范围为20′(约36 km)、30′(约54 km)、40′(约72 km)的4°×6°范围区域地形改正进行计算，待求点个数均为3 110 400个，结果见表 3。

表 3 不同积分半径严格棱柱体积分串行及并行计算耗时Tab. 3 Serial and parallel computing of strict prism method with different integral radius

积分半径/(′)	积分半径内待积分点数	严格棱柱体积分串行计算耗时/h	严格棱柱积分八分量法并行计算耗时/min
20	58 080	49.927	3.766
30	130 320	114.221	8.422
40	231 360	208.266	14.633

表选项

由表 3可以看出，当地形分辨率较高时(分辨率为10″)，随着积分半径的增加，棱柱法计算局部地形改正的耗时呈线性增长。利用传统棱柱法积分计算不同积分半径范围的耗时均在2 d以上，本文所提八分量并行算法则可有效提高计算效率，计算时间均在15 min以内，充分利用了计算设备的计算能力，使得在有限时间内计算高分辨率局部地形改正成为可能。

3  全国范围地形改正快速计算与结果分析

通过上述并行策略与算法改进，使得在短时间内计算高分辨率局部地形改正成为可能。结合上述方法，本节对中国陆地主要区域(南海周边岛屿未完全涵盖)局部地形改正进行计算，地形格网分辨率为30″，积分半径选取2°，计算范围由18°N至54°N，跨度36°；由72°E至136°E，跨度64°。并行算法计算耗时199.47 min，串行计算耗时超过70 d。计算结果如表 4、图 6所示。

表 4 中国陆地主要区域地形改正数值特征Tab. 4 Numerical characteristics of topographic correction in main land areas of ChinamGal

最大值	最小值	平均值	中位数
196.30	0.00	1.831	0.257

表选项

图 6 中国陆地主要区域地形改正 Fig. 6     Topographic correction map of major land areas in China      

图选项

由表 4、图 6可知，本文所提并行算法可以有效解决传统棱柱法对于大范围高分辨率地形改正计算耗时过长的缺点，在200 min以内即可完成对中国陆地主要区域30″的地形改正计算。由计算结果可知我国地形改正量普遍低于80 mGal，平均改正值1.83 mGal，改正值中位数0.257 mGal。改正量较大值出现在喜马拉雅山脉、横断山等地形起伏剧烈地区，改正量最大值达到196 mGal。西部地区由于山脉较多，地形起伏相对剧烈，改正量较大，平均改正值高于30 mGal。东部平原地区由于地势平坦，改正量普遍低于10 mGal。台湾地区由于中央山脉等高山影响，地形改正量相对较大，最高达到61.36 mGal。

4  结论与讨论

本文利用CPU+GPU的异构并行方法可以较好地实现大范围区域地形改正快速并行计算。经过试验计算分析表明：

(1) 本文所提严格棱柱八分量并行算法，可以很好地适应GPU计算规则与设备运算逻辑。通过地形格网重新编码和严格棱柱体积分拆解方法，有效解决了计算量分配与线程计算能力间的矛盾，实现了基于GPU+CPU异构并行技术的严格棱柱体积分计算地形改正快速并行算法。在利用CPU+GPU异构并行计算时，本文方法可以充分利用GPU设备性能，加速比最高达到850倍以上，有效提高了计算效率。进行4°×6°范围，积分半径40′，分辨率1′的局部地形改正计算仅需1.5 s；分辨率10″的局部地形改正计算仅需14.6 min；进行分辨率3″地形改正计算耗时45.7 h，传统串行算法则难以完成计算。

(2) 精度方面。本文算法与传统严格棱柱体积分串行算法间的最大绝对误差可保证在10^-6 mGal以内，平均误差和均方根误差均小于10^-12 mGal，计算精度较高。而近区使用线性卷积分的快速算法和严格棱柱体积分串行算法相比，最大误差达到4.77 mGal，平均差值为0.11 mGal；RMS值为0.32 mGal。因此本文算法无论是效率还是精度均优于传统串行快速算法。

(3) 利用本文算法对我国陆地主要区域进行地形改正计算，结果表明：我国地形改正量值普遍低于80 mGal，平均改正值1.83 mGal。改正量较大值出现在喜马拉雅山脉、横断山等地形起伏剧烈地区，改正量最大值达到196 mGal。西部地区由于山脉较多，地形起伏相对剧烈，改正量较大，平均改正值高于30 mGal；东部平原地区由于地势平坦，改正量普遍低于10 mGal；台湾地区由于中央山脉等高山影响，地形改正量相对较大，最高达到61.36 mGal。

作者简介

第一作者简介：黄炎(1994-), 男, 博士生, 研究方向为物理大地测量。E-mail:[e-mail]

通信作者：冯进凯, E-mail：[e-mail]

权威 | 专业 | 学术 | 前沿

微信、抖音小视频投稿邮箱 | [e-mail]

微信公众号中搜索「测绘学报」，关注我们，长按上图二维码，关注学术前沿动态。

欢迎加入《测绘学报》作者QQ群： 751717395

进群请备注：姓名+单位+稿件编号